6 休息中の凝集性に影響する要因
以下の分析では、どのような要因が群れ内の凝集性に影響しているかを検討します。
6.1 データの加工
まず、休息中の周辺個体数などについて算出して元データに結合します。
RG_female: 休息集団内のメス数
RG_female_plus1: 自身を含めた休息集団内のメス数
RG_est: 休息集団内の発情メス数
RG_**: それぞれのオスが休息集団内にいたか
x3m_female: 3m以内のメス数
x5m_female: 5m以内のメス数
**_3m: それぞれのオスが3m以内にいるか
## 各調査期間ごとの6歳以上のメスのID
adult18 <- unique(female_18m %>% filter(age >= 6) %>% .$femaleID)
adult19 <- unique(female_19m %>% filter(age >= 6) %>% .$femaleID)
adult20 <- unique(female_20m %>% filter(age >= 6) %>% .$femaleID)
adult21 <- unique(female_21m %>% filter(age >= 6) %>% .$femaleID)
## 10m以内近接個個体名を記したデータを作成
focal_prox <- focal_combined_final %>%
separate(x0_1m, into = str_c("x0_1m",1:11), sep = ",") %>%
separate(x1_3m, into = str_c("x1_3m",1:15), sep = ",") %>%
separate(x3_5m, into = str_c("x3_5m",1:16), sep = ",") %>%
separate(x5_10m, into = str_c("x5_10m", 1:9), sep = ",") %>%
pivot_longer(cols = x0_1m1:x5_10m9,
names_to = "proximity",
values_to = "ID") %>%
mutate(proximity = ifelse(str_detect(proximity,"x0_1m"),"x0_1m",
ifelse(str_detect(proximity,"x1_3m"),"x1_3m",
ifelse(str_detect(proximity,"x3_5m"),"x3_5m",
ifelse(str_detect(proximity,"x5_10m"),"x5_10m","NA"))))) %>%
filter(!is.na(ID))
## 休息集団に含まれるメスの数
RG_female <- focal_prox %>%
filter(RG == "1") %>%
filter((study_period == "m18" & ID %in% adult18)|(study_period == "nm19" & ID %in% adult18)|
(study_period == "m19" & ID %in% adult19)|
(study_period == "m20" & ID %in% adult20)|(study_period == "nm21" & ID %in% adult20)|
(study_period == "m21" & ID %in% adult21)|(study_period == "nm22" & ID %in% adult21)) %>%
left_join(female_all %>% select(date, femaleID, rs2),
by = c("date","ID" = "femaleID")) %>%
group_by(date, no_focal, time) %>%
summarise(RG_female = n(),
RG_est = sum(rs2, na.rm = TRUE)) %>%
ungroup()
## 休息集団内に各オスがいるか否か
RG_TY <- focal_prox %>%
filter(RG == "1") %>%
select(date, no_focal, time, ID) %>%
filter(ID == "TY") %>%
mutate(RG_TY = 1) %>%
select(-ID)
RG_IT <- focal_prox %>%
filter(RG == "1") %>%
select(date, no_focal, time, ID) %>%
filter(ID == "IT") %>%
mutate(RG_IT = 1) %>%
select(-ID)
RG_LK <- focal_prox %>%
filter(RG == "1") %>%
select(date, no_focal, time, ID) %>%
filter(ID == "LK") %>%
mutate(RG_LK = 1) %>%
select(-ID)
RG_KR <- focal_prox %>%
filter(RG == "1") %>%
select(date, no_focal, time, ID) %>%
filter(ID == "KR") %>%
mutate(RG_KR = 1) %>%
select(-ID)
## 3m、5m近接(相手個体の活動に依らず)
focal_prox_all <- focal_prox %>%
mutate(ID = str_replace(ID, "\\(F\\)","")) %>%
mutate(ID = str_replace(ID, "\\(M\\)","")) %>%
mutate(ID = str_replace(ID, "\\(O\\)","")) %>%
mutate(ID = str_replace(ID, "\\(N\\)","")) %>%
mutate(ID = str_replace(ID, "\\(n\\)","")) %>%
mutate(ID = str_replace(ID, "\\(LT\\)","")) %>%
filter((study_period == "m18" & ID %in% adult18)|(study_period == "nm19" & ID %in% adult18)|
(study_period == "m19" & ID %in% adult19)|
(study_period == "m20" & ID %in% adult20)|(study_period == "nm21" & ID %in% adult20)|
(study_period == "m21" & ID %in% adult21)|(study_period == "nm22" & ID %in% adult21))
x3m_female <- focal_prox_all %>%
filter(proximity %in% c("x0_1m","x1_3m")) %>%
group_by(date, no_focal, time) %>%
summarise(x3m_female = n()) %>%
ungroup()
x5m_female <- focal_prox_all %>%
filter(proximity %in% c("x0_1m","x1_3m","x3_5m")) %>%
group_by(date, no_focal, time) %>%
summarise(x5m_female = n()) %>%
ungroup()
## 元データに結合
focal_combined_prox <- focal_combined_final %>%
left_join(RG_female, by = c("date","no_focal","time")) %>%
left_join(x3m_female, by = c("date","no_focal","time")) %>%
left_join(x5m_female, by = c("date","no_focal","time")) %>%
left_join(RG_TY, by = c("date","no_focal","time")) %>%
left_join(RG_IT, by = c("date","no_focal","time")) %>%
left_join(RG_LK, by = c("date","no_focal","time")) %>%
left_join(RG_KR, by = c("date","no_focal","time")) %>%
replace_na(list(RG_female = 0, RG_est = 0, x3m_female = 0, x5m_female =0, RG_TY = 0,
RG_IT = 0, RG_LK = 0, RG_KR = 0)) %>%
replace_na(list(x0_1m = "NA",x1_3m = "NA", x3_5m = "NA", x5_10m = "NA")) %>%
## 自身を含んだ休息集団サイズ
mutate(RG_female_plus1 = RG_female + 1) %>%
## オスとの近接情報を追加
mutate(TY_3m = ifelse(str_detect(x0_1m,"TY")|str_detect(x1_3m,"TY"),1,0),
IT_3m = ifelse(str_detect(x0_1m,"IT")|str_detect(x1_3m,"IT"),1,0),
KR_3m = ifelse(str_detect(x0_1m,"KR")|str_detect(x1_3m,"KR"),1,0),
LK_3m = ifelse(str_detect(x0_1m,"LK")|str_detect(x1_3m,"LK"),1,0)) 続いて、その日群れ内にいた血縁個体数を計算し、結合します。
kin <- read_csv("../Data/data/others/kin.csv")
female_presence <- group_all %>%
select(groupID, study_period, date, Kur:Yun) %>%
select(-c(TY,IT, LK, KR, KM, TG)) %>%
pivot_longer(cols = Kur:Yun,
names_to = "femaleID",
values_to = "presence") %>%
left_join(att) %>%
## 6歳以上の個体のみを抽出
filter(age >= 6,
presence == 1)
focal_numkin <- focal_combined_final %>%
distinct(study_period, date, groupID, no_focal, subject) %>%
left_join(female_presence %>% select(date, groupID, femaleID),
by = c("date", "groupID")) %>%
filter(subject != femaleID) %>%
left_join(kin, by = c("subject" = "femaleID", "femaleID" = "femaleID2")) %>%
mutate(kin = ifelse(kin >= 0.0625,1,0)) %>%
group_by(no_focal, date, subject) %>%
summarise(no_kin = sum(kin)) %>%
ungroup()
focal_combined_prox %>%
left_join(focal_numkin, by = c("no_focal","date","subject")) -> focal_combined_prox_b続いて、TY、ITとの親密度(CSI)の情報も結合します。
focal_combined_prox_b %>%
left_join(CSI_TY_combined %>% select(subject, CSI_TY)) %>%
left_join(CSI_IT_combined %>% select(subject, CSI_IT)) -> focal_combined_prox_c最後に、その日観察したオトナメスへの攻撃頻度も追加します。
### 6歳以上への攻撃のみ
aggression_all %>%
mutate(femaleID = str_replace_all(femaleID, "\\?","")) %>%
left_join(att, by = c("femaleID", "study_period")) %>%
## 確実に被攻撃個体が6歳以上であるものを抽出
filter(age >= 6 | str_detect(femaleID,"multi|many|Multi|KunTrt|Ako,Kil")) %>%
group_by(date) %>%
summarise(no_agg = length(unique(no))) %>%
ungroup() %>%
filter(date >= "2018-10-08") -> no_agg_6yo
focal_combined_prox_c %>%
left_join(group_all_b %>%
select(groupID, date, duration)) %>%
left_join(no_agg_6yo, by = "date") %>%
mutate(rate_agg_6yo = no_agg*60/duration) %>%
replace_na(list(rate_agg_6yo = 0)) %>%
left_join(female_all %>% select(date, femaleID, rs2),
by = c("subject" = "femaleID",
"date")) -> focal_combined_prox_final作成したデータは以下の通り。
6.2 休息集団サイズの分析
ここでは、TYやITの有無を含む様々な要因が個体追跡中の休息集団サイズ(正確には、休息集団内のメス数)に影響するかを検討します。
6.2.1 交尾期(非発情メス)
6.2.1.1 データの加工
まず、非発情メスについて分析を行うためのデータを作成します。連続変数は標準化します。
- 休息ポイント数が10以上の個体追跡セッションのみを使用
- ハドリングのポイントは除去
- 地上休息のみを含む
cohesion_data_an <- focal_combined_prox_final %>%
mutate(rate_female = no_female/max_female) %>%
filter(rs2 == 0) %>%
replace_na(list(hud = 0)) %>%
filter(hud != "1") %>%
filter(RG == 1 & T_G == "G") %>%
group_by(date, no_focal, subject, study_period, rate_agg_6yo, no_ntm,
no_est, no_male, no_female, TY_presence,
no_kin, IT_presence, CSI_TY) %>%
summarise(sum_RG = sum(RG_female_plus1),
n = n()) %>%
ungroup() %>%
mutate(TY_presence = as.factor(TY_presence),
IT_presence = as.factor(IT_presence)) %>%
filter(study_period != "m18") %>%
filter(n >= 10) %>%
mutate(ntm_std = standardize(no_ntm),
male_std = standardize(no_male),
est_std = standardize(no_est),
kin_std = standardize(no_kin),
female_std = standardize(no_female),
rate_agg_std = standardize(rate_agg_6yo),
CSI_TY_std = standardize(CSI_TY),
mean_RG = sum_RG/n)
cohesion_data_an共変量と応答変数(平均休息集団サイズ)の関連を図示すると以下のようになります。
cohesion_data_an %>%
pivot_longer(cols = c(no_ntm, no_est, no_female, rate_agg_6yo, kin_std, CSI_TY),
names_to = "type",
values_to = "value") %>%
ggplot(aes(x = value, y = mean_RG)) +
geom_point(shape = 1) +
facet_rep_wrap(~ type,
scales = "free") +
theme_bw(base_size = 12) +
theme(aspect.ratio = 1) +
labs(x = "", y = "Mean resting group size")
6.2.1.2 モデリング
それでは、モデリングを行います。以下のように行いました。連続変数は標準化しています。また、ITの効果を見る際には、2021年交尾期を除くモデルを作成しました。
- 応答変数: log(平均休息集団サイズ(メスのみ))
- 分布: Studentのt分布
- 説明変数: TYの在不在(
TY_presence)、ITの在不在(IT_presence)、群れ外オス数(ntm_std)、集団内の発情メス数(est_std)、集団内のメス数(female_std)、集団内の血縁メス数(kin_est)、調査期間(study_period)
- ランダム切片:
subject
## 2019~2021年
m_cohesion_an <- brm(log(mean_RG) ~ TY_presence + IT_presence + ntm_std + est_std +
female_std + kin_std + study_period + (1|subject),
family = student,
iter = 11000, warmup = 1000, seed = 123,
prior = c(prior(student_t(4,0,5), class = "b"),
prior(student_t(4,0,10), class = "Intercept"),
prior(student_t(4,0,10), class = "sd"),
prior(student_t(4,0,10), class = "nu"),
prior(student_t(4,0,10), class = "sigma")),
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "cmdstanr",
data = cohesion_data_an,
file = "model/m_cohesion_an.rds")
## 2019~2020年
cohesion_data_an_bf21 <- cohesion_data_an %>%
filter(study_period != "m21") %>%
mutate(ntm_std = standardize(no_ntm),
male_std = standardize(no_male),
est_std = standardize(no_est),
kin_std = standardize(no_kin),
female_std = standardize(no_female),
rate_agg_std = standardize(rate_agg_6yo),
CSI_TY_std = standardize(CSI_TY),
mean_RG = sum_RG/n)
m_cohesion_an_bf21 <- brm(log(mean_RG) ~ TY_presence + IT_presence + ntm_std + est_std +
female_std + kin_std + study_period + (1|subject),
family = student,
iter = 11000, warmup = 1000, seed = 123,
prior = c(prior(student_t(4,0,5), class = "b"),
prior(student_t(4,0,10), class = "Intercept"),
prior(student_t(4,0,10), class = "sd")),
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "cmdstanr",
data = cohesion_data_an_bf21,
file = "model/m_cohesion_an_bf21.rds")6.2.1.3 モデルチェック
まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。
多重共線性の問題もありません。
bayesplotパッケージ(Gabry & Mahr, 2022)のpp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。
Rhatにも問題はなく、収束の問題はないよう。
data.frame(Rhat = brms::rhat(m_cohesion_an)) %>%
ggplot(aes(x = Rhat))+
geom_histogram(fill = "white",
color = "black")+
theme_bw()+
theme(aspect.ratio = 1)
data.frame(Rhat = brms::rhat(m_cohesion_an_bf21)) %>%
ggplot(aes(x = Rhat))+
geom_histogram(fill = "white",
color = "black")+
theme_bw()+
theme(aspect.ratio = 1)
6.2.1.4 結果の確認
まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。6.2.1.5 結果の図示
最後に、結果の図示をします。
まず、TYの有無についてです。
emm_cohesion_an_TY <- emmeans(m_cohesion_an,
~TY_presence)
cohesion_data_an %>%
ggplot(aes(x = TY_presence, y = log(mean_RG))) +
geom_violin(bw = 0.2,
scale = "width") +
geom_boxplot(size = 0.5,
linewidth = 0.4,
width = 0.2,
fill = "grey",
outlier.alpha = 0) +
geom_point(data = emm_cohesion_an_TY %>%
data.frame() %>%
mutate(TY_presence = as.factor(TY_presence)),
aes(y = emmean),
size = 3,
color = "white") +
geom_errorbar(data = emm_cohesion_an_TY %>%
data.frame() %>%
mutate(TY_presence = as.factor(TY_presence)),
aes(y = emmean, ymin = lower.HPD,
ymax = upper.HPD),
size = 0.4,
color = "white",
width = 0.1) +
geom_signif(comparisons = list(c("0", "1")),
y_position = 2.5,
annotations = "***") +
theme_bw(base_size = 12) +
scale_x_discrete(labels = c("*TY* absent", "*TY* present")) +
labs(x = "", y = "log(mean resting group size)") +
theme(text = element_text(family = "Times New Roman"),
axis.text.x = element_markdown(),
aspect.ratio = 1) -> p_cohesion_an_TY
p_cohesion_an_TY
ggsave("figures/p_cohesion_an_TY.tif", p_cohesion_an_TY,
dpi = 2000, units = "mm",
width = 100, height = 100)続いて、ITについてです。
emm_cohesion_an_IT <- emmeans(m_cohesion_an_bf21,
~IT_presence)
cohesion_data_an_bf21 %>%
ggplot(aes(x = IT_presence, y = log(mean_RG))) +
geom_violin(bw = 0.2,
scale = "width") +
geom_boxplot(size = 0.5,
linewidth = 0.4,
width = 0.2,
fill = "grey",
outlier.alpha = 0) +
geom_point(data = emm_cohesion_an_IT %>%
data.frame() %>%
mutate(IT_presence = as.factor(IT_presence)),
aes(y = emmean),
size = 3,
color = "white") +
geom_errorbar(data = emm_cohesion_an_IT %>%
data.frame() %>%
mutate(IT_presence = as.factor(IT_presence)),
aes(y = emmean, ymin = lower.HPD,
ymax = upper.HPD),
size = 0.4,
color = "white",
width = 0.1) +
theme_bw(base_size = 12) +
scale_x_discrete(labels = c("*IT* absent", "*IT* present")) +
labs(x = "", y = "log(mean resting group size)") +
theme(text = element_text(family = "Times New Roman"),
axis.text.x = element_markdown(),
aspect.ratio = 1) -> p_cohesion_an_IT
p_cohesion_an_IT
ggsave("figures/p_cohesion_an_IT.tif", p_cohesion_an_IT,
dpi = 2000, units = "mm",
width = 100, height = 100)6.2.2 非交尾期
6.2.2.1 データの加工
続いて、非交尾期について分析を行うためのデータを作成します。連続変数は標準化します。
- 休息ポイント数が10以上の個体追跡セッションのみを使用
- ハドリングのポイントは除去
- 地上休息のみを含む
cohesion_data_nm <- focal_combined_prox_final %>%
mutate(rate_female = no_female/max_female) %>%
filter(str_detect(study_period, "nm")) %>%
replace_na(list(hud = 0)) %>%
filter(hud != "1") %>%
filter(RG == 1 & T_G == "G") %>%
group_by(date, no_focal, subject, study_period,
no_female, TY_presence,
no_kin, IT_presence) %>%
summarise(sum_RG = sum(RG_female_plus1),
n = n()) %>%
ungroup() %>%
mutate(TY_presence = as.factor(TY_presence),
IT_presence = as.factor(IT_presence)) %>%
filter(study_period != "m18") %>%
filter(n >= 10) %>%
mutate(kin_std = standardize(no_kin),
female_std = standardize(no_female),
mean_RG = sum_RG/n)
cohesion_data_nm共変量と応答変数(平均休息集団サイズ)の関連を図示すると以下のようになります。
cohesion_data_nm %>%
pivot_longer(cols = c(no_female, no_kin),
names_to = "type",
values_to = "value") %>%
ggplot(aes(x = value, y = mean_RG)) +
geom_point(shape = 1) +
facet_rep_wrap(~ type,
scales = "free") +
theme_bw(base_size = 12) +
theme(aspect.ratio = 1) +
labs(x = "", y = "Mean resting group size")
6.2.2.2 モデリング
それでは、モデリングを行います。以下のように行いました。連続変数は標準化しています。また、ITの効果を見る際には、2021年交尾期を除くモデルを作成しました。
- 応答変数: log(平均休息集団サイズ(メスのみ))
- 分布: Studentのt分布
- 説明変数: TYの在不在(
TY_presence)、ITの在不在(IT_presence)、群れ外オス数(ntm_std)、集団内の発情メス数(est_std)、集団内のメス数(female_std)、集団内の血縁メス数(kin_est)、調査期間(study_period)
- ランダム切片:
subject
## 2019~2021年
m_cohesion_nm <- brm(log(mean_RG) ~ TY_presence + kin_std + study_period + (1|subject),
family = student,
iter = 11000, warmup = 1000, seed = 123,
prior = c(prior(student_t(4,0,5), class = "b"),
prior(student_t(4,0,10), class = "Intercept"),
prior(student_t(4,0,10), class = "sd"),
prior(student_t(4,0,10), class = "nu"),
prior(student_t(4,0,10), class = "sigma")),
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "cmdstanr",
data = cohesion_data_nm,
file = "model/m_cohesion_nm.rds")
## 切断モデル
m_cohesion_nm_trunc <- brm(log(mean_RG)|trunc(lb = 0) ~ TY_presence +
kin_std + study_period + (1|subject),
family = student,
iter = 7000, warmup = 2000, seed = 123,
prior = c(prior(student_t(4,0,5), class = "b"),
prior(student_t(4,0,10), class = "Intercept"),
prior(student_t(4,0,10), class = "sd"),
prior(student_t(4,0,10), class = "nu"),
prior(student_t(4,0,10), class = "sigma")),
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "cmdstanr",
data = cohesion_data_nm,
file = "model/m_cohesion_nm.rds")
## 2019~2020年
cohesion_data_nm_bf21 <- cohesion_data_nm %>%
filter(study_period != "m22") %>%
mutate(kin_std = standardize(no_kin),
female_std = standardize(no_female),
mean_RG = sum_RG/n)
m_cohesion_nm_bf21 <- brm(log(mean_RG) ~ TY_presence + kin_std + study_period + (1|subject),
family = student,
iter = 11000, warmup = 1000, seed = 123,
prior = c(prior(student_t(4,0,5), class = "b"),
prior(student_t(4,0,10), class = "Intercept"),
prior(student_t(4,0,10), class = "sd")),
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "cmdstanr",
data = cohesion_data_nm_bf21,
file = "model/m_cohesion_nm_bf21.rds")
## 切断モデル
m_cohesion_nm_bf21_trunc <- brm(log(mean_RG)|trunc(lb = 0) ~ TY_presence + IT_presence + ntm_std + est_std +
female_std + kin_std + study_period + (1|subject),
family = student,
iter = 11000, warmup = 1000, seed = 123,
prior = c(prior(student_t(4,0,5), class = "b"),
prior(student_t(4,0,10), class = "Intercept"),
prior(student_t(4,0,10), class = "sd")),
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "cmdstanr",
data = cohesion_data_nm_bf21,
file = "model/m_cohesion_nm_bf21.rds")6.2.2.3 モデルチェック
まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。
多重共線性の問題もありません。
bayesplotパッケージ(Gabry & Mahr, 2022)のpp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。
Rhatにも問題はなく、収束の問題はないよう。
data.frame(Rhat = brms::rhat(m_cohesion_nm)) %>%
ggplot(aes(x = Rhat))+
geom_histogram(fill = "white",
color = "black")+
theme_bw()+
theme(aspect.ratio = 1)
data.frame(Rhat = brms::rhat(m_cohesion_nm_bf21)) %>%
ggplot(aes(x = Rhat))+
geom_histogram(fill = "white",
color = "black")+
theme_bw()+
theme(aspect.ratio = 1)
6.2.2.4 結果の確認
まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。6.2.2.5 結果の図示
最後に、結果の図示をします。
まず、TYの有無についてです。
emm_cohesion_nm_TY <- emmeans(m_cohesion_nm,
~TY_presence)
cohesion_data_nm %>%
ggplot(aes(x = TY_presence, y = log(mean_RG))) +
geom_violin(bw = 0.2,
scale = "width") +
geom_boxplot(size = 0.5,
linewidth = 0.4,
width = 0.2,
fill = "grey",
outlier.alpha = 0) +
geom_point(data = emm_cohesion_nm_TY %>%
data.frame() %>%
mutate(TY_presence = as.factor(TY_presence)),
aes(y = emmean),
size = 3,
color = "white") +
geom_errorbar(data = emm_cohesion_nm_TY %>%
data.frame() %>%
mutate(TY_presence = as.factor(TY_presence)),
aes(y = emmean, ymin = lower.HPD,
ymax = upper.HPD),
size = 0.4,
color = "white",
width = 0.1) +
theme_bw(base_size = 12) +
scale_x_discrete(labels = c("*TY* absent", "*TY* present")) +
labs(x = "", y = "log(mean resting group size)") +
theme(text = element_text(family = "Times New Roman"),
axis.text.x = element_markdown(),
aspect.ratio = 1) -> p_cohesion_nm_TY
p_cohesion_nm_TY
ggsave("figures/p_cohesion_nm_TY.tif", p_cohesion_nm_TY,
dpi = 2000, units = "mm",
width = 100, height = 100)6.3 近接個体数の分析(休息中)
ここでは、TYやITの有無を含む様々な要因が個体追跡中の周辺メス数に影響するかを検討します。
6.3.1 交尾期(非発情メス)
6.3.1.1 データの加工
まず、非発情メスについて分析を行うためのデータを作成します。連続変数は標準化します。
- 休息ポイント数が10以上の個体追跡セッションのみを使用
- ハドリングのポイントは除去
- 地上休息のみを含む
prox_data_an_RG <- focal_combined_prox_final %>%
mutate(rate_female = no_female/max_female) %>%
filter(rs2 == 0) %>%
replace_na(list(hud = 0)) %>%
filter(hud != "1") %>%
filter(RG == 1 & T_G == "G") %>%
group_by(date, no_focal, subject, study_period, rate_agg_6yo, no_ntm,
no_est, no_male, no_female, TY_presence,
no_kin, IT_presence, CSI_TY) %>%
summarise(sum_prox = sum(x3m_female + 1),
n = n()) %>%
ungroup() %>%
mutate(TY_presence = as.factor(TY_presence),
IT_presence = as.factor(IT_presence)) %>%
filter(study_period != "m18") %>%
filter(n >= 10) %>%
mutate(ntm_std = standardize(no_ntm),
male_std = standardize(no_male),
est_std = standardize(no_est),
kin_std = standardize(no_kin),
female_std = standardize(no_female),
rate_agg_std = standardize(rate_agg_6yo),
CSI_TY_std = standardize(CSI_TY),
mean_prox = sum_prox/n)
prox_data_an_RG共変量と応答変数(平均休息集団サイズ)の関連を図示すると以下のようになります。
prox_data_an_RG %>%
pivot_longer(cols = c(no_ntm, no_est, no_female, rate_agg_6yo, kin_std, CSI_TY),
names_to = "type",
values_to = "value") %>%
ggplot(aes(x = value, y = mean_prox)) +
geom_point(shape = 1) +
facet_rep_wrap(~ type,
scales = "free") +
theme_bw(base_size = 12) +
theme(aspect.ratio = 1) +
labs(x = "", y = "Mean resting group size")
6.3.1.2 モデリング
それでは、モデリングを行います。以下のように行いました。連続変数は標準化しています。また、ITの効果を見る際には、2021年交尾期を除くモデルを作成しました。
- 応答変数: log(平均休息集団サイズ(メスのみ))
- 分布: Studentのt分布
- 説明変数: TYの在不在(
TY_presence)、ITの在不在(IT_presence)、群れ外オス数(ntm_std)、集団内の発情メス数(est_std)、集団内のメス数(female_std)、集団内の血縁メス数(kin_est)、調査期間(study_period)
- ランダム切片:
subject
## 2019~2021年
m_prox_an_RG <- brm(log(mean_prox) ~ TY_presence + IT_presence + ntm_std + est_std +
female_std + kin_std + study_period + (1|subject) + (1|date),
family = student,
iter = 11000, warmup = 1000, seed = 123,
prior = c(prior(student_t(4,0,5), class = "b"),
prior(student_t(4,0,10), class = "Intercept"),
prior(student_t(4,0,10), class = "sd"),
prior(student_t(4,0,10), class = "nu"),
prior(student_t(4,0,10), class = "sigma")),
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "cmdstanr",
data = prox_data_an_RG,
file = "model/m_prox_an_RG.rds")
## 2019~2020年
prox_data_an_RG_bf21 <- prox_data_an_RG %>%
filter(study_period != "m21") %>%
mutate(ntm_std = standardize(no_ntm),
male_std = standardize(no_male),
est_std = standardize(no_est),
kin_std = standardize(no_kin),
female_std = standardize(no_female),
rate_agg_std = standardize(rate_agg_6yo),
CSI_TY_std = standardize(CSI_TY),
mean_prox = sum_prox/n)
m_prox_an_RG_bf21 <- brm(log(mean_prox) ~ TY_presence + IT_presence + ntm_std + est_std +
female_std + kin_std + study_period + (1|subject),
family = student,
iter = 11000, warmup = 1000, seed = 123,
prior = c(prior(student_t(4,0,5), class = "b"),
prior(student_t(4,0,10), class = "Intercept"),
prior(student_t(4,0,10), class = "sd")),
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "cmdstanr",
data = prox_data_an_RG_bf21,
file = "model/m_prox_an_RG_bf21.rds")6.3.1.3 モデルチェック
まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。
多重共線性の問題もありません。
bayesplotパッケージ(Gabry & Mahr, 2022)のpp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。
Rhatにも問題はなく、収束の問題はないよう。
data.frame(Rhat = brms::rhat(m_prox_an_RG)) %>%
ggplot(aes(x = Rhat))+
geom_histogram(fill = "white",
color = "black")+
theme_bw()+
theme(aspect.ratio = 1)
data.frame(Rhat = brms::rhat(m_prox_an_RG_bf21)) %>%
ggplot(aes(x = Rhat))+
geom_histogram(fill = "white",
color = "black")+
theme_bw()+
theme(aspect.ratio = 1)
6.3.1.4 結果の確認
まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。6.3.1.5 結果の図示
最後に、結果の図示をします。
まず、TYの有無についてです。
emm_prox_an_RG_TY <- emmeans(m_prox_an_RG,
~TY_presence)
prox_data_an_RG %>%
ggplot(aes(x = TY_presence, y = log(mean_prox))) +
geom_violin(bw = 0.2,
scale = "width") +
geom_boxplot(size = 0.5,
linewidth = 0.4,
width = 0.2,
fill = "grey",
outlier.alpha = 0) +
geom_point(data = emm_prox_an_RG_TY %>%
data.frame() %>%
mutate(TY_presence = as.factor(TY_presence)),
aes(y = emmean),
size = 3,
color = "white") +
geom_errorbar(data = emm_prox_an_RG_TY %>%
data.frame() %>%
mutate(TY_presence = as.factor(TY_presence)),
aes(y = emmean, ymin = lower.HPD,
ymax = upper.HPD),
size = 0.4,
color = "white",
width = 0.1) +
geom_signif(comparisons = list(c("0", "1")),
y_position = 2.1,
annotations = "***") +
theme_bw(base_size = 12) +
scale_x_discrete(labels = c("*TY* absent", "*TY* present")) +
labs(x = "", y = "log(mean no. of females within 3m)") +
theme(text = element_text(family = "Times New Roman"),
axis.text.x = element_markdown(),
aspect.ratio = 1) -> p_prox_an_RG_TY
p_prox_an_RG_TY
ggsave("figures/p_prox_an_RG_TY.tif", p_prox_an_RG_TY,
dpi = 2000, units = "mm",
width = 100, height = 100)続いて、ITについてです。
emm_prox_an_RG_IT <- emmeans(m_prox_an_RG_bf21,
~IT_presence)
prox_data_an_RG_bf21 %>%
ggplot(aes(x = IT_presence, y = log(mean_prox))) +
geom_violin(bw = 0.2,
scale = "width") +
geom_boxplot(size = 0.5,
linewidth = 0.4,
width = 0.2,
fill = "grey",
outlier.alpha = 0) +
geom_point(data = emm_prox_an_RG_IT %>%
data.frame() %>%
mutate(IT_presence = as.factor(IT_presence)),
aes(y = emmean),
size = 3,
color = "white") +
geom_errorbar(data = emm_prox_an_RG_IT %>%
data.frame() %>%
mutate(IT_presence = as.factor(IT_presence)),
aes(y = emmean, ymin = lower.HPD,
ymax = upper.HPD),
size = 0.4,
color = "white",
width = 0.1) +
theme_bw(base_size = 12) +
scale_x_discrete(labels = c("*IT* absent", "*IT* present")) +
labs(x = "", y = "log(mean resting group size)") +
theme(text = element_text(family = "Times New Roman"),
axis.text.x = element_markdown(),
aspect.ratio = 1) -> p_prox_an_RG_IT
p_prox_an_RG_IT
ggsave("figures/p_prox_an_RG_IT.tif", p_prox_an_RG_IT,
dpi = 2000, units = "mm",
width = 100, height = 100)6.3.2 非交尾期
6.3.2.1 データの加工
続いて、非交尾期について分析を行うためのデータを作成します。連続変数は標準化します。
- 休息ポイント数が10以上の個体追跡セッションのみを使用
- ハドリングのポイントは除去
- 地上休息のみを含む
prox_data_nm_RG <- focal_combined_prox_final %>%
mutate(rate_female = no_female/max_female) %>%
filter(str_detect(study_period, "nm")) %>%
replace_na(list(hud = 0)) %>%
filter(hud != "1") %>%
filter(RG == 1 & T_G == "G") %>%
group_by(date, no_focal, subject, study_period,
no_female, TY_presence,
no_kin, IT_presence) %>%
summarise(sum_prox = sum(x3m_female + 1),
n = n()) %>%
ungroup() %>%
mutate(TY_presence = as.factor(TY_presence),
IT_presence = as.factor(IT_presence)) %>%
filter(study_period != "m18") %>%
filter(n >= 10) %>%
mutate(kin_std = standardize(no_kin),
female_std = standardize(no_female),
mean_prox = sum_prox/n)
prox_data_nm_RG共変量と応答変数(平均休息集団サイズ)の関連を図示すると以下のようになります。
prox_data_nm_RG %>%
pivot_longer(cols = c(no_female, no_kin),
names_to = "type",
values_to = "value") %>%
ggplot(aes(x = value, y = mean_prox)) +
geom_point(shape = 1) +
facet_rep_wrap(~ type,
scales = "free") +
theme_bw(base_size = 12) +
theme(aspect.ratio = 1) +
labs(x = "", y = "Mean resting group size")
6.3.2.2 モデリング
それでは、モデリングを行います。以下のように行いました。連続変数は標準化しています。また、ITの効果を見る際には、2021年交尾期を除くモデルを作成しました。
- 応答変数: log(平均休息集団サイズ(メスのみ))
- 分布: Studentのt分布
- 説明変数: TYの在不在(
TY_presence)、ITの在不在(IT_presence)、群れ外オス数(ntm_std)、集団内の発情メス数(est_std)、集団内のメス数(female_std)、集団内の血縁メス数(kin_est)、調査期間(study_period)
- ランダム切片:
subject
## 2019~2021年
m_prox_nm_RG <- brm(log(mean_prox) ~ TY_presence + kin_std + study_period + (1|subject),
family = student,
iter = 11000, warmup = 1000, seed = 123,
prior = c(prior(student_t(4,0,5), class = "b"),
prior(student_t(4,0,10), class = "Intercept"),
prior(student_t(4,0,10), class = "sd"),
prior(student_t(4,0,10), class = "nu"),
prior(student_t(4,0,10), class = "sigma")),
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "cmdstanr",
data = prox_data_nm_RG,
file = "model/m_prox_nm_RG.rds")6.3.2.3 モデルチェック
まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。
多重共線性の問題もありません。
bayesplotパッケージ(Gabry & Mahr, 2022)のpp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。
Rhatにも問題はなく、収束の問題はないよう。
data.frame(Rhat = brms::rhat(m_prox_nm_RG)) %>%
ggplot(aes(x = Rhat))+
geom_histogram(fill = "white",
color = "black")+
theme_bw()+
theme(aspect.ratio = 1)
6.3.2.5 結果の図示
最後に、結果の図示をします。
まず、TYの有無についてです。
emm_prox_nm_RG_TY <- emmeans(m_prox_nm_RG,
~TY_presence)
prox_data_nm_RG %>%
ggplot(aes(x = TY_presence, y = log(mean_prox))) +
geom_violin(bw = 0.2,
scale = "width") +
geom_boxplot(size = 0.5,
linewidth = 0.4,
width = 0.2,
fill = "grey",
outlier.alpha = 0) +
geom_point(data = emm_prox_nm_RG_TY %>%
data.frame() %>%
mutate(TY_presence = as.factor(TY_presence)),
aes(y = emmean),
size = 3,
color = "white") +
geom_errorbar(data = emm_prox_nm_RG_TY %>%
data.frame() %>%
mutate(TY_presence = as.factor(TY_presence)),
aes(y = emmean, ymin = lower.HPD,
ymax = upper.HPD),
size = 0.4,
color = "white",
width = 0.1) +
theme_bw(base_size = 12) +
scale_x_discrete(labels = c("*TY* absent", "*TY* present")) +
labs(x = "", y = "log(mean resting group size)") +
theme(text = element_text(family = "Times New Roman"),
axis.text.x = element_markdown(),
aspect.ratio = 1) -> p_prox_nm_RG_TY
p_prox_nm_RG_TY
ggsave("figures/p_prox_nm_RG_TY.tif", p_prox_nm_RG_TY,
dpi = 2000, units = "mm",
width = 100, height = 100)6.4 近接個体数の分析(採食中)
最後に、TYやITの有無を含む様々な要因が個体追跡中の周辺メス数に影響するかを検討します。
6.4.1 交尾期(非発情メス)
6.4.1.1 データの加工
まず、非発情メスについて分析を行うためのデータを作成します。連続変数は標準化します。
- 休息ポイント数が10以上の個体追跡セッションのみを使用
- ハドリングのポイントは除去
- 地上休息のみを含む
prox_data_an_F <- focal_combined_prox_final %>%
mutate(rate_female = no_female/max_female) %>%
filter(rs2 == 0) %>%
replace_na(list(hud = 0)) %>%
filter(hud != "1") %>%
filter(activity == "F" & T_G == "G") %>%
group_by(date, no_focal, subject, study_period, rate_agg_6yo, no_ntm,
no_est, no_male, no_female, TY_presence,
no_kin, IT_presence, CSI_TY) %>%
summarise(sum_prox = sum(x3m_female + 1),
n = n()) %>%
ungroup() %>%
mutate(TY_presence = as.factor(TY_presence),
IT_presence = as.factor(IT_presence)) %>%
filter(study_period != "m18") %>%
filter(n >= 10) %>%
mutate(ntm_std = standardize(no_ntm),
male_std = standardize(no_male),
est_std = standardize(no_est),
kin_std = standardize(no_kin),
female_std = standardize(no_female),
rate_agg_std = standardize(rate_agg_6yo),
CSI_TY_std = standardize(CSI_TY),
mean_prox = sum_prox/n)
prox_data_an_F共変量と応答変数(平均休息集団サイズ)の関連を図示すると以下のようになります。
prox_data_an_F %>%
pivot_longer(cols = c(no_ntm, no_est, no_female, rate_agg_6yo, kin_std, CSI_TY),
names_to = "type",
values_to = "value") %>%
ggplot(aes(x = value, y = mean_prox)) +
geom_point(shape = 1) +
facet_rep_wrap(~ type,
scales = "free") +
theme_bw(base_size = 12) +
theme(aspect.ratio = 1) +
labs(x = "", y = "Mean resting group size")
6.4.1.2 モデリング
それでは、モデリングを行います。以下のように行いました。連続変数は標準化しています。また、ITの効果を見る際には、2021年交尾期を除くモデルを作成しました。
- 応答変数: log(平均休息集団サイズ(メスのみ))
- 分布: Studentのt分布
- 説明変数: TYの在不在(
TY_presence)、ITの在不在(IT_presence)、群れ外オス数(ntm_std)、集団内の発情メス数(est_std)、集団内のメス数(female_std)、集団内の血縁メス数(kin_est)、調査期間(study_period)
- ランダム切片:
subject
## 2019~2021年
m_prox_an_F <- brm(log(mean_prox) ~ TY_presence + IT_presence + ntm_std + est_std +
female_std + kin_std + study_period + (1|subject) + (1|date),
family = student,
iter = 11000, warmup = 1000, seed = 123,
prior = c(prior(student_t(4,0,5), class = "b"),
prior(student_t(4,0,10), class = "Intercept"),
prior(student_t(4,0,10), class = "sd"),
prior(student_t(4,0,10), class = "nu"),
prior(student_t(4,0,10), class = "sigma")),
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "cmdstanr",
data = prox_data_an_F,
file = "model/m_prox_an_F.rds")
## 2019~2020年
prox_data_an_F_bf21 <- prox_data_an_F %>%
filter(study_period != "m21") %>%
mutate(ntm_std = standardize(no_ntm),
male_std = standardize(no_male),
est_std = standardize(no_est),
kin_std = standardize(no_kin),
female_std = standardize(no_female),
rate_agg_std = standardize(rate_agg_6yo),
CSI_TY_std = standardize(CSI_TY),
mean_prox = sum_prox/n)
m_prox_an_F_bf21 <- brm(log(mean_prox) ~ TY_presence + IT_presence + ntm_std + est_std +
female_std + kin_std + study_period + (1|subject),
family = student,
iter = 11000, warmup = 1000, seed = 123,
prior = c(prior(student_t(4,0,5), class = "b"),
prior(student_t(4,0,10), class = "Intercept"),
prior(student_t(4,0,10), class = "sd")),
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "cmdstanr",
data = prox_data_an_F_bf21,
file = "model/m_prox_an_F_bf21.rds")6.4.1.3 モデルチェック
まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。
多重共線性の問題もありません。
bayesplotパッケージ(Gabry & Mahr, 2022)のpp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。
Rhatにも問題はなく、収束の問題はないよう。
data.frame(Rhat = brms::rhat(m_prox_an_F)) %>%
ggplot(aes(x = Rhat))+
geom_histogram(fill = "white",
color = "black")+
theme_bw()+
theme(aspect.ratio = 1)
data.frame(Rhat = brms::rhat(m_prox_an_F_bf21)) %>%
ggplot(aes(x = Rhat))+
geom_histogram(fill = "white",
color = "black")+
theme_bw()+
theme(aspect.ratio = 1)
6.4.1.4 結果の確認
まずは、全期間のデータを用いたモデルについてみる。太字になっている変数は95%確信区間が0をまたいでおらず、有意な影響があるといえる。6.4.1.5 結果の図示
最後に、結果の図示をします。
まず、TYの有無についてです。
emm_prox_an_F_TY <- emmeans(m_prox_an_F,
~TY_presence)
prox_data_an_F %>%
ggplot(aes(x = TY_presence, y = log(mean_prox))) +
geom_violin(bw = 0.1,
scale = "width") +
geom_boxplot(size = 0.5,
linewidth = 0.4,
width = 0.2,
fill = "grey",
outlier.alpha = 0) +
geom_point(data = emm_prox_an_F_TY %>%
data.frame() %>%
mutate(TY_presence = as.factor(TY_presence)),
aes(y = emmean),
size = 3,
color = "white") +
geom_errorbar(data = emm_prox_an_F_TY %>%
data.frame() %>%
mutate(TY_presence = as.factor(TY_presence)),
aes(y = emmean, ymin = lower.HPD,
ymax = upper.HPD),
size = 0.4,
color = "white",
width = 0.1) +
geom_signif(comparisons = list(c("0", "1")),
y_position = 1.3,
annotations = "***") +
theme_bw(base_size = 12) +
scale_x_discrete(labels = c("*TY* absent", "*TY* present")) +
labs(x = "", y = "log(mean no. of females within 3m)") +
theme(text = element_text(family = "Times New Roman"),
axis.text.x = element_markdown(),
aspect.ratio = 1) -> p_prox_an_F_TY
p_prox_an_F_TY
ggsave("figures/p_prox_an_F_TY.tif", p_prox_an_F_TY,
dpi = 2000, units = "mm",
width = 100, height = 100)続いて、ITについてです。
emm_prox_an_F_IT <- emmeans(m_prox_an_F_bf21,
~IT_presence)
prox_data_an_F_bf21 %>%
ggplot(aes(x = IT_presence, y = log(mean_prox))) +
geom_violin(bw = 0.07,
scale = "width") +
geom_boxplot(size = 0.5,
linewidth = 0.4,
width = 0.2,
fill = "grey",
outlier.alpha = 0) +
geom_point(data = emm_prox_an_F_IT %>%
data.frame() %>%
mutate(IT_presence = as.factor(IT_presence)),
aes(y = emmean),
size = 3,
color = "white") +
geom_errorbar(data = emm_prox_an_F_IT %>%
data.frame() %>%
mutate(IT_presence = as.factor(IT_presence)),
aes(y = emmean, ymin = lower.HPD,
ymax = upper.HPD),
size = 0.4,
color = "white",
width = 0.1) +
theme_bw(base_size = 12) +
scale_x_discrete(labels = c("*IT* absent", "*IT* present")) +
labs(x = "", y = "log(mean resting group size)") +
theme(text = element_text(family = "Times New Roman"),
axis.text.x = element_markdown(),
aspect.ratio = 1) -> p_prox_an_F_IT
p_prox_an_F_IT
ggsave("figures/p_prox_an_F_IT.tif", p_prox_an_F_IT,
dpi = 2000, units = "mm",
width = 100, height = 100)6.4.2 非交尾期
6.4.2.1 データの加工
続いて、非交尾期について分析を行うためのデータを作成します。連続変数は標準化します。
- 休息ポイント数が10以上の個体追跡セッションのみを使用
- ハドリングのポイントは除去
- 地上休息のみを含む
prox_data_nm_F <- focal_combined_prox_final %>%
mutate(rate_female = no_female/max_female) %>%
filter(str_detect(study_period, "nm")) %>%
replace_na(list(hud = 0)) %>%
filter(hud != "1") %>%
filter(activity == "F" & T_G == "G") %>%
group_by(date, no_focal, subject, study_period,
no_female, TY_presence,
no_kin, IT_presence) %>%
summarise(sum_prox = sum(x3m_female + 1),
n = n()) %>%
ungroup() %>%
mutate(TY_presence = as.factor(TY_presence),
IT_presence = as.factor(IT_presence)) %>%
filter(study_period != "m18") %>%
filter(n >= 10) %>%
mutate(kin_std = standardize(no_kin),
female_std = standardize(no_female),
mean_prox = sum_prox/n)
prox_data_nm_F共変量と応答変数(平均休息集団サイズ)の関連を図示すると以下のようになります。
prox_data_nm_F %>%
pivot_longer(cols = c(no_female, no_kin),
names_to = "type",
values_to = "value") %>%
ggplot(aes(x = value, y = mean_prox)) +
geom_point(shape = 1) +
facet_rep_wrap(~ type,
scales = "free") +
theme_bw(base_size = 12) +
theme(aspect.ratio = 1) +
labs(x = "", y = "Mean resting group size")
6.4.2.2 モデリング
それでは、モデリングを行います。以下のように行いました。連続変数は標準化しています。また、ITの効果を見る際には、2021年交尾期を除くモデルを作成しました。
- 応答変数: log(平均休息集団サイズ(メスのみ))
- 分布: Studentのt分布
- 説明変数: TYの在不在(
TY_presence)、ITの在不在(IT_presence)、群れ外オス数(ntm_std)、集団内の発情メス数(est_std)、集団内のメス数(female_std)、集団内の血縁メス数(kin_est)、調査期間(study_period)
- ランダム切片:
subject
## 2019~2021年
m_prox_nm_F <- brm(log(mean_prox) ~ TY_presence + kin_std + study_period + (1|subject),
family = student,
iter = 11000, warmup = 1000, seed = 195,
prior = c(prior(student_t(4,0,5), class = "b"),
prior(student_t(4,0,10), class = "Intercept"),
prior(student_t(4,0,10), class = "sd"),
prior(student_t(4,0,10), class = "nu"),
prior(student_t(4,0,10), class = "sigma")),
control = list(adapt_delta = 0.99, max_treedepth = 20),
backend = "cmdstanr",
data = prox_data_nm_F,
file = "model/m_prox_nm_F.rds")6.4.2.3 モデルチェック
まず、DHARMaパッケージ(Hartig, 2022)とDHARMa.helperパッケージ(Rodríguez-Sánchez, 2023)でモデルの前提が満たされているかを確認します。
多重共線性の問題もありません。
bayesplotパッケージ(Gabry & Mahr, 2022)のpp_check関数で、事後分布からの予測分布と実測値の分布を比較しても大きな乖離はない。
Rhatにも問題はなく、収束の問題はないよう。
data.frame(Rhat = brms::rhat(m_prox_nm_F)) %>%
ggplot(aes(x = Rhat))+
geom_histogram(fill = "white",
color = "black")+
theme_bw()+
theme(aspect.ratio = 1)
6.4.2.5 結果の図示
最後に、結果の図示をします。
まず、TYの有無についてです。
emm_prox_nm_F_TY <- emmeans(m_prox_nm_F,
~TY_presence)
prox_data_nm_F %>%
ggplot(aes(x = TY_presence, y = log(mean_prox))) +
geom_violin(bw = 0.2,
scale = "width") +
geom_boxplot(size = 0.5,
linewidth = 0.4,
width = 0.2,
fill = "grey",
outlier.alpha = 0) +
geom_point(data = emm_prox_nm_F_TY %>%
data.frame() %>%
mutate(TY_presence = as.factor(TY_presence)),
aes(y = emmean),
size = 3,
color = "white") +
geom_errorbar(data = emm_prox_nm_F_TY %>%
data.frame() %>%
mutate(TY_presence = as.factor(TY_presence)),
aes(y = emmean, ymin = lower.HPD,
ymax = upper.HPD),
size = 0.4,
color = "white",
width = 0.1) +
theme_bw(base_size = 12) +
scale_x_discrete(labels = c("*TY* absent", "*TY* present")) +
labs(x = "", y = "log(mean resting group size)") +
theme(text = element_text(family = "Times New Roman"),
axis.text.x = element_markdown(),
aspect.ratio = 1) -> p_prox_nm_F_TY
p_prox_nm_F_TY
ggsave("figures/p_prox_nm_F_TY.tif", p_prox_nm_F_TY,
dpi = 2000, units = "mm",
width = 100, height = 100)References
Gabry, J., & Mahr, T. (2022). Bayesplot: Plotting for bayesian models. https://mc-stan.org/bayesplot
Hartig, F. (2022). DHARMa: Residual diagnostics for hierarchical (multi-level / mixed) regression models. https://CRAN.R-project.org/package=DHARMa
Rodríguez-Sánchez, F. (2023). DHARMa.helpers: Helper functions to check models not (yet) directly supported by DHARMa.